/*
 *  MWRK.cpp
 *  T3nsors
 *
 *  Created by Michael Barriault on 10-06-10.
 *  Copyright 2010 University of Guelph. All rights reserved.
 *
 */

#include "MW.h"


void MW::RK2(double dt) {
	Scalar Phi2,Psi2;
	Vector E2,B2;
	Phi2 = Psi2 = Scalar(O);
	E2 = B2 = Vector(O,3);
	
	Scalar k1,k2;
	Scalar l1,l2;
	Vector m1,m2;
	Vector n1,n2;
	
	k1 = Phidot(Phi,Psi,E,B);
	l1 = Psidot(Phi,Psi,E,B);
	m1 = Edot(Phi,Psi,E,B);
	n1 = Bdot(Phi,Psi,E,B);
	
	Boundaries(&k1,&l1,&m1,&n1);
	
	FORO Phi2[o] = Phi[o] + dt*k1[o];
	FORO Psi2[o] = Psi[o] + dt*l1[o];
	FORa FORO E2(a)[o] = E(a)[o] + dt*m1(a)[o];
	FORa FORO B2(a)[o] = B(a)[o] + dt*n1(a)[o];
	
	k2 = Phidot(Phi2,Psi2,E2,B2);
	l2 = Psidot(Phi2,Psi2,E2,B2);
	m2 = Edot(Phi2,Psi2,E2,B2);
	n2 = Bdot(Phi2,Psi2,E2,B2);
	
	Boundaries(&k2,&l2,&m2,&n2);
	
	FORO Phi[o] = Phi[o] + dt*(k1[o]+k2[o])/2;
	FORO Psi[o] = Psi[o] + dt*(l1[o]+l2[o])/2;
	FORa FORO E(a)[o] = E(a)[o] + dt*(m1(a)[o]+m2(a)[o])/2;
	FORa FORO B(a)[o] = B(a)[o] + dt*(n1(a)[o]+n2(a)[o])/2;
	
	return;
}

void MW::RK3(double dt) {
	Scalar Phi2,Phi3,Psi2,Psi3;
	Vector E2,E3,B2,B3;
	Phi2 = Phi3 = Psi2 = Psi3 = Scalar(O);
	E2 = E3 = B2 = B3 = Vector(O,3);
	
	Scalar k1,k2,k3;
	Scalar l1,l2,l3;
	Vector m1,m2,m3;
	Vector n1,n2,n3;
	
	k1 = Phidot(Phi,Psi,E,B);
	l1 = Psidot(Phi,Psi,E,B);
	m1 = Edot(Phi,Psi,E,B);
	n1 = Bdot(Phi,Psi,E,B);
	
	Boundaries(&k1,&l1,&m1,&n1);
	
	FORO Phi2[o] = Phi[o] + 0.5*dt*k1[o];
	FORO Psi2[o] = Psi[o] + 0.5*dt*l1[o];
	FORa FORO E2(a)[o] = E(a)[o] + 0.5*dt*m1(a)[o];
	FORa FORO B2(a)[o] = B(a)[o] + 0.5*dt*n1(a)[o];
	
	k2 = Phidot(Phi2,Psi2,E2,B2);
	l2 = Psidot(Phi2,Psi2,E2,B2);
	m2 = Edot(Phi2,Psi2,E2,B2);
	n2 = Bdot(Phi2,Psi2,E2,B2);
	
	Boundaries(&k2,&l2,&m2,&n2);
	
	FORO Phi3[o] = Phi[o] + 0.75*dt*k2[o];
	FORO Psi3[o] = Psi[o] + 0.75*dt*l2[o];
	FORa FORO E3(a)[o] = E(a)[o] + 0.75*dt*m2(a)[o];
	FORa FORO B3(a)[o] = B(a)[o] + 0.75*dt*n2(a)[o];
	
	k3 = Phidot(Phi3,Psi3,E3,B3);
	l3 = Psidot(Phi3,Psi3,E3,B3);
	m3 = Edot(Phi3,Psi3,E3,B3);
	n3 = Bdot(Phi3,Psi3,E3,B3);
	
	Boundaries(&k3,&l3,&m3,&n3);
	
	FORO Phi[o] = Phi[o] + dt*(2*k1[o]+3*k2[o]+4*k3[o])/9;
	FORO Psi[o] = Psi[o] + dt*(2*l1[o]+3*l2[o]+4*l3[o])/9;
	FORa FORO E(a)[o] = E(a)[o] + dt*(2*m1(a)[o]+3*m2(a)[o]+4*m3(a)[o])/9;
	FORa FORO B(a)[o] = B(a)[o] + dt*(2*n1(a)[o]+3*n2(a)[o]+4*n3(a)[o])/9;
	
	return;
}

